The following starts on restructuring the linear predictor into a multi-level form. The revision effects are constructed from a overall, joint and assignment level effects. The backbone duration will follow a similar template.
Below we consider parameters that characterise variation in the linear predictor due to:
Running MCMC with 2 parallel chains...
Chain 2 finished in 18.6 seconds.
Chain 1 finished in 20.0 seconds.
Both chains finished successfully.
Mean chain execution time: 19.3 seconds.
Total execution time: 20.0 seconds.
Warning: 12 of 2000 (1.0%) transitions ended with a divergence.
See https://mc-stan.org/misc/warnings for details.
Fit model to simulated data
# quick work - but estimates will be rough, variance very rough# snk <- capture.output(# f1 <- m1$pathfinder(ld, num_paths=20, single_path_draws=200,# history_size=50, max_lbfgs_iters=100,# refresh = 0, draws = 2000)# )
Comparing parameter estimates with parameters used to simulate the data.
Estimating the log-odds of response averaged across the sample, comparing observed vs model.
Log-odds response (surgery)
d_post <-data.table(f1$draws(variables =c("eta"), format ="matrix"))d_post <-melt(d_post, measure.vars =names(d_post))d_post[, i :=gsub("eta[","",variable,fixed=T)]d_post[, i :=as.numeric(gsub("]","",i,fixed=T))]d_fig <-copy(d)d_fig[, i :=1:.N]d_post <-merge(d_post, d_fig, by ="i")# create a new variable to indicate what was receivedd_post[, s_d1 :=copy(u_d1)]# set all that assigned to dair to have recvd daird_post[s ==2& d1 ==0, s_d1 :=1]d_fig <- d_post[ , .(mu =mean(value), eta =mean(eta)), keyby = .(s,j,d1,s_d1)]d_fig[, s :=factor(s, labels =c("early", "late", "chronic"))]d_fig[, j :=factor(j, labels =c("hip", "knee"))]d_fig[, s_d1 :=factor(s_d1, labels =c("dair", "one", "two"))]d_fig[, d1 :=factor(d1, labels =c("not-rand", "rand-dair", "rand-rev"))]ggplot(d_fig, aes(x = d1, y = mu, col = s_d1)) +geom_point(aes(x = d1, y = eta, col = s_d1), pch =3, size=2) +geom_point(size =1) +scale_x_discrete("Assigned/selected trt arm") +scale_y_continuous("log-odds response") +scale_color_discrete("Selected surgery") +facet_grid(j ~ s, scales ="free_x")
Figure 2: Estimated log-odds of treatment success vs true
Decomposition of revision effects into group means, within group means and effects. Prior and posterior view on group variation - what you can take from this is that the variance components are poorly informed by the data, i.e. there isn’t going to be any pooling because there are not enough groups to learn from.
Variance components (revision)
d_fig <-data.table( f1$draws(variables =c("s_mu_d1", "s_b_d1"), format ="matrix"))d_fig <-melt(d_fig, measure.vars =names(d_fig))ggplot(d_fig, aes(x = value, group = variable)) +geom_density() +stat_function(fun = fGarch::dstd, args =list(nu =3, mean =0, sd =2), col =2, lty =2) +facet_wrap(~variable)
Figure 3: Between and within SD
Revision effects
# effectsd_fig <-data.table( f1$draws(variables =c("b_d1[5,1]", "b_d1[6,1]","b_d1[5,2]", "b_d1[6,2]"), format ="matrix"))d_fig <-melt(d_fig, measure.vars =names(d_fig))d_fig <- d_fig[, .(mu_b_d1 =mean(value)), keyby = variable]# overall meand_fig <-cbind( d_fig, data.table( f1$draws(variables =c("mu_d1"), format ="matrix"))[, .(mu_d1_rev =mean(mu_d1))] )# joint level meand_post <-data.table(f1$draws(variables =c("mu_d1_j"), format ="matrix"))d_fig[, mu_d1_rev_jnt :=rep(colMeans(d_post), each =2)]# jointd_fig[, j :=gsub("b_d1\\[[5-6],", "", variable)]d_fig[, j :=gsub("\\],", "", j)]d_fig[, j :=factor(j, labels =c("hip", "knee"))]# rev type (one/two)d_fig[, type :=gsub("b_d1\\[", "", variable)]d_fig[, type :=substr(type, 1, 1)]d_fig[, type :=factor(type, labels =c("one", "two"))]# on any given small sample, probably isn't going to be that close.# d_fig[, tru_mu_d1_rev := par$mu_d1_rev]# d_fig[j == "hip", tru_mu_d1_rev_j := par$mu_d1_rev_jnt[1]]# d_fig[j == "knee", tru_mu_d1_rev_j := par$mu_d1_rev_jnt[2]]ggplot(d_fig, aes(x = variable, y = mu_b_d1)) +geom_point() +geom_hline(aes(yintercept = mu_d1_rev)) +geom_hline(aes(yintercept = mu_d1_rev_jnt), lty =2) +# geom_hline(aes(yintercept = tru_mu_d1_rev), col = 2) +# geom_hline(aes(yintercept = tru_mu_d1_rev_j), lty = 2, col = 2) +scale_x_discrete("Assigned/selected trt arm") +scale_y_continuous("log-odds ratio") +facet_wrap(j + type ~., nrow =1, scales ="free_x")
Figure 4: Revision effects, between, within and assignment level means